APS/123-QED 

Constrained Monte Carlo Method and 
Calculation of the Temperature Dependence of Magnetic 

Anisotropy 

P. Asselin 

Seagate Technology, Bloomington, MN 55435, USA 

R. F. L. Evans0 J. Barker, and R. W. Chantrell 

Department of Physics, University of York, 
Heslington, York YOlO 5DD United Kingdom 

R. Yanes and O. Chubykalo-Fesenko 
Instituto de Ciencia de Materiales de Madrid, 
CSIC, Cantoblanco, Madrid 28049, Spain 

D. Hinzke and U. Nowak 
Universitat Konstanz, Fachbereich Physik, 
Universitdtsstrafie 10, D-78464 Konstanz, Germany 
(Dated: June 18, 2010) 



1 



Abstract 

We introduce a constrained Monte Carlo method which allows us to traverse the phase space 
of a classical spin system while fixing the magnetization direction. Subsequently we show the 
method's capability to model the temperature dependence of magnetic anisotropy, and for bulk 
uniaxial and cubic anisotropics we recover the low-temperature Callen-Callen power laws in M. 
We also calculate the temperature scaling of the 2-ion anisotropy in LIq FePt, and recover the 
experimentally observed M^-^ scaling. The method is newly applied to evaluate the temperature- 
dependent effective anisotropy in the presence of the Neel surface anisotropy in thin films with 
different easy axis configurations. In systems having different surface and bulk easy axes, we show 
the capability to model the temperature-induced reorientation transition. The intrinsic surface 
anisotropy is found to follow a linear temperature behavior in a large range of temperatures. 

PACS numbers: 75.30.Gw, 75.70.Rf, 75.10.Hk, 75.70.Ak 
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I. INTRODUCTION 



The temperature dependence of magnetocrystalline anisotropy of pure ferromagnets has 
been well known for decades, following the Callen-Callen theory-. Examples of such systems 
classed as pure ferromagnets in the anisotropic sense include Gadolinium (Uniaxial) and Fe 
(Cubic). Other ferromagnets, such as Co, have a much more complicated temperature 
dependence, due to the crystallographic origin of the anisotropy^. Indeed, with increased 
temperature the anisotropy in Co exhibits a change in sign, indicating a transformation 
from an easy-axis to easy-plane anisotropy. Such behavior is not explained by the Callen- 
Callen theory. Other materials not exhibiting a simple temperature dependence of the 
anisotropy are magnetic transition metal alloys such as FePt and CoPt. Here the origin of the 
anisotropy is due to an ion-ion anisotropic exchange interaction, arising from the underlying 
crystal symmetry. In FePt the anisotropy exhibits an unusual temperature dependence^i^ 
of oc M^'^. Other systems of technological interest include magnetic thin films and 
nanoparticles, where surface effects can lead to unusual temperature dependent anisotropics. 
In general the temperature dependence of anisotropy in many materials is not obvious, and 
as such is still an active area of research some 40 years after the Callen-Callen theory. 

Recently, the high temperature behavior of magnetic anisotropy has become important 
due to the applications in heat-assisted magnetic recording (HAMR)^""-. The idea of HAMR 
is based on the heating of the recording media to decrease the writing field of the high 
anisotropy media (such as FePt) to values compatible with the writing fields provided by 
conventional recording heads. Since the writing field is proportional to the anisotropy field 
iJk = 2K^{T)/M{T), the knowledge of the scaling behavior of the anisotropy with the 
magnetization M has become a paramount consideration for HAMR^. It should be noted 
that even in relatively simple systems, a simple scaling behavior predicted by the Callen- 
Callen theory is only valid at temperatures far from the Curie temperature. The systems 
proposed for HAMR applications can also include more complex composite media such as 
soft/hard bilayers^, FePt/FeRh with metamagnetic phase transition^°''^^, or exchange-bias 
systemai^. 

The evaluation of the temperature dependence of magnetic anisotropy is also important 
for the modeling of the laser-induced demagnetization processes. The thermal decrease of 
the anisotropy during the laser-induced demagnetization has been shown to be responsi- 
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ble for the optically-induced magnetization precession^^. Thus the ability to evaluate the 
temperature dependence of the anisotropy in complex systems at arbitrary temperatures is 
highly desired from the fundamental and applied perspectives. 

In this sense magnetic thin films with surface anisotropy are a representative example of 
this more complicated situation. In ultra-thin films, especially when in contact with a dif- 
ferent non-magnetic matrix, the interplay between the broken symmetry, magnetostriction, 
roughness, spin-orbit interaction and charge transfer can often be encompassed in a phe- 
nomenological model as an additional surface anisotropy. Since the surface anisotropy has 
a different temperature dependence from the bulk, multiple experiments on thin film have 
demonstrated the occurrence of the spin reorientation transition from an out-of-plane to 
in-plane magnetization both as a function of temperature and thin film thickness^"—. The 
possibility to engineer the reorientation transition also requires the capability to evaluate 
the temperature dependence of the surface anisotropy independently from the bulk. 

In the following we present a new Monte Carlo method which can be applied to the com- 
putation of both bulk and surface anisotropics at finite temperature. This article represents 
a first step showing the possibility to calculate the temperature-dependent anisotropics in 
principle. When combined with detailed magnetic information, such as that available from 
ab-initio methods^- , this forms a very powerful method of engineering the temperature de- 
pendent properties of a magnetic system. 

II. MODELING METHODS 

For the calculations presented in the following we describe the magnetic properties of 
the system by utilizing a classical atomistic spin model, similar to Nowak— , with a general 
Hamiltonian of the form: 



describing the exchange interaction {Jij), uniaxial (-ft'u), cubic (Kc), and Neel surface 
anisotropics^^ {Ks) respectively. Note that in the following text, we also refer to the ef- 
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fective temperature dependent values (uniaxial), Kf^ (cubic), and K^^ (surface). The 
summations for the exchange and surface anisotropy are generally limited to nearest neigh- 
bors only, except for the case of FePt where the full exchange up to 5 neighboring cells 
(around 1300 neighbors) was taken into account. 

The parameters are chosen to represent a generic ferromagnet, with a Curie 
temperature (Tc) of around lOOOK, and arbitrary anisotropy constants, where K^i,Kc <^ Jij. 
Note that all three anisotropy terms have been included within the same Hamiltonian for 
brevity - in practice a system will only have cubic or uniaxial anisotropy, and clearly only 
surface atoms will possess surface anisotropy. 

We compute thermodynamic properties by averaging over the Boltzmann distribution 
using the Metropolis algorithnt^. Our innovation, which we call the constrained Monte Carlo 
method (CMC), is to modify the elementary moves of the random walk so as to conserve 
the average magnetization direction M = (Xli ^^^^ sample the 

Boltzmann distribution over a submanifold of the full phase space. Thus we keep the system 
out of thermodynamic equilibrium in a controlled manner, while allowing its microscopic 
degrees of freedom to thermalize. 

Because the system cannot reach full equilibrium, the average of the total internal torque 
T = ^ &H/dSi') does not vanish. We show in Appendix [Cl that this is equal to 

the macroscopic torque — M x 9J-'/9M, where J^(M) is the Helmholtz free energy, now a 
function of M. Even though we cannot compute J-" directly, we can reconstruct its angular 
dependence by integration. 



where the integral on M' can be taken along any path on which the system behaves 
reversibly. This in turn gives us the anisotropy constants at any temperature. 

In practice it is often simpler to recover the anisotropy constants directly from the deriva- 
tives. We first initialize the system with uniform magnetization in a direction of our choice, 
away from the anisotropy axes, where we expect a nonzero torque. Next we evolve the sys- 
tem by constrained Monte Carlo until the length of the magnetization reaches equilibrium. 
We then take a thermodynamic average of the torque over a large number constrained Monte 
Carlo steps, typically 50,000. We repeat at other orientations and we finally reconstruct the 
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anisotropy constants from the angular dependence of the torque. 

III. CONSTRAINED MONTE CARLO 

The Metropohs algorithm works by generating trial moves at random and accepting or 
rejecting each move based on the ratio of the Boltzmann probability densities exp(— 
P = 1/kT, at the initial and final states. This ratio depends only on the energy difference 
between the two states. An accepted move yields a new state and a rejected move yields a 
repetition of the initial state ("null move" in the following). There is considerable freedom in 
the construction of trial moves. It is required only that each move have the same probability 
density as the inverse move (reversibility), and that all states be reachable by a sequence of 
moves (ergodicity). Under these conditions the random walk's limiting distribution is the 
Boltzmann distribution. 

In the usual Monte Carlo method we generate our trial moves by drawing a vector v 
from an isotropic normal distribution, choosing a spin Sj at random, adding Sj x v to it and 
normalizing the result to obtain a trial spin S^. The probability density of the move depends 
only on the angle between Sj and S^, which ensures reversibility. Ergodicity is obvious. The 
variance of the v distribution controls the size of the attempted moves and can be chosen 
at will to improve the ratio of accepted to rejected moves, similarly to the parameter a in 
reference |26|. For our Hamiltonian, the energy difference involves only spin i and the few 
neighboring spins to which it is coupled by exchange, so the decision to accept or reject the 
move can be made quickly. A sequence of moves, counting null moves, constitutes a step; 
we compute quantities of interest once per step to average them. 

In the constrained Monte Carlo method the trial moves act on two spins at a time. The 
extra degrees of freedom allow us to fix M to any given unit vector, which we take here to 
be the positive z axis since we can always reduce the problem to this case by means of a 
global rotation. Ignoring the z coordinates for a moment, we simply displace two spins by 
equal and opposite amounts in the XY plane. There are technical issues due to the fact 
are not canonical coordinates; furthermore, we need to allow sign changes in the 
z coordinates. In the end we settled on the following: 

1. Choose a primary spin Sj and a compensation spin Sj, not necessarily neighbors. 
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2. Displace the primary spin as in the usual Monte Carlo method, obtaining a new spin 

3. Adjust the compensation spin's x and y components to preserve = My = 0, 



q' — q. _i_ Q. — q' 

q' — q. A- q. _ q' 

jy jy ^ ^y iy 



4. Adjust the z component, 



s'^, = sign^s,,)^i-s';i-s';i . 

If the argument of the square root is negative, stop and take a null move. 

5. Compute the new magnetization. 

If M'^ < 0, stop and take a null move. 

6. Compute the energy difference ATi = 71' — H. 

7. Compute the acceptance probability P, 



P = min 



8. Accept the move with probability P or take a null move with probability 1 — P. 

In effect, we use the compensation spin to project the system back to its admissible 
manifold. The projection is not orthogonal and does not preserve measure. Consequently 
the Boltzmann ratio in step [7] is multiplied by a geometric correction, the ratio of two 
Jacobians, which we derive in Appendix |X1 We prove ergodicity in Appendix |Bl 

The null moves at step H] handle the kinematic constraints in a natural way. We could 
instead add fictitious states that allow steplHto complete and assign zero probability (infinite 
energy) to these states; this would guarantee rejection at step [8] and it is simpler to stop at 
step m Similarly at step |5] we reject trial move that would change the sign of M; here the 
end states exist but we simply want to sample them with zero probability. 
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IV. CALCULATION OF BULK ANISOTROPIES 



Given the originality of the constrained Monte Carlo method, it is important to ensure 
that the method is reliable and conforms to existing results, especially regarding the low 
temperature dependence of bulk anisotropy as predicted by Callen and Calleni, where uni- 
axial anisotropy was shown to have an dependence, and cubic anisotropy to have an 
dependence. 

In the present work such bulk systems were approximated by simulating a generic ferro- 
magnetic system with 16000 spin moments with periodic boundary conditions, to eliminate 
surface effects and minimize finite size effects. Sample torque curves for the system with 
uniaxial anisotropy are shown in Fig. [TJ The points show the calculated torque and the 
curves are fits to a sin (26') dependence, where 6 is the angle from the easy axis. The sin (2^) 
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FIG. 1: Simulated and analytical angular dependence of the restoring torque for uniaxial 
anisotropy at temperatures 10 K, 100 K, and 500 K. (Color Online) 



relationship is seen to hold at all temperatures and the fitted proportionality constant gives 
K^(T). In a situation such as this, where all the torque curves have the same shape, the 
anisotropy is described by a single parameter and it is sufficient to compute the torque at 
the maximum. 

In more general cases, such as those with surface anisotropy, it is necessary to compute 
the torque at several angular positions. Finally, with every new system it is prudent to 
verify the shape of the torque curves over many angles, both polar and azimuthal, before 
reducing the number of points to the minimum necessary. 
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The uniaxial anisotropy, cubic anisotropy and magnetization for the generic system are 



plotted against temperature in Fig. 2(a) In order to reduce the computational effort, the 



torques were computed at a single angular position, 9 = 45° in the uniaxial case and 6 = 22.5° 



in the x2;-plane for cubic anisotropy. In Fig. 2(b) the anisotropies are plotted against the 



magnetization on logarithmic scales. As can be seen, the results show excellent agreement 
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(a)Temperature dependence of magnetization, uniaxial, 
and cubic anisotropies. The lines provide a guide to the 

eye. 
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(b)Temperature scaling of uniaxial and cubic 
anisotropies. 

FIG. 2: Simulation results for the temperature dependence and scaling of a pure ferromagnet with 
uniaxial and cubic anisotropies. (Color Online) 



at low-temperatures with the scaling behavior predicted by Callen and Callen. We note that 
in these calculations the internal energy and free energy are almost interchangeable. That 
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is, dT /do is nearly equal to d{H)/d6. The anisotropy is too weak to affect the entropy and 
the difference (H) — T = TS is nearly independent of 6. We expect this to hold in most 
magnetic systems, but in full generality the distinction between J-" and ("H) must be kept. 

In strong anisotropy systems and temperatures close to we have found that the system 
torque deviates from the expected sin (2^) behavior. This is in agreement with previous 
publications where it has been shown that at high temperatures strong magnetization fluc- 
tuations lead to several interesting effects, related to the free energy behavior, such as the 
change of the magnetization length on the saddle poiniM or the elliptic character of domain 

V. CALCULATION OF ANISOTROPY IN FePt 

FePt is a material of current interest because of its extremely high magnetocrystalline 
anisotropy energy in the LIq crystal phase^^. The material is also unusual due to the M^'^ 
low-temperature scaling of the effective anisotropy^*^. Ab-initio and Langevin dynamics 
simulations by Mryasov et alM managed to reproduce the observed scaling of the anisotropy 
by calculation of the internal energy, and in this work we have applied the Constrained Monte 
Carlo method to the same problem. We have utilized the same Hamiltonian as Mryasov et 
alM and simulated a system 6 nm^ with periodic boundary conditions. The key addition to 
the generic Hamiltonian in Eq. ([T]) is a two-ion anisotropy term of the form: 



where K^j is the anisotropy constant which is site and range dependent, as extracted form 
the ab-initio calculations. The system also possesses an easy-plane anisotropy which is 
an order of magnitude weaker than the easy-axis 2-ion anisotropy. The existence of these 
competing anisotropics in fact gives rise to the unusual scaling exponent due to the different 
temperature scaling of the anisotropic contributions, for two-ion and for single ion. 
Due to the large value of the effective anisotropy in FePt, the torque deviates from the 
expected sin (2^) relationship at elevated temperatures, and so the free energy was obtained 
by integration over 6. The torque curves were calculated by first equilibrating the system for 
10000 steps at each temperature and angle, and then the thermal average of the torque was 
calculated over a further 70000 steps. Plots of the change in free energy, AF, are shown in 




(3) 
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FIG. 3: Angular dependence of the free energy in LIq FePt at high temperatures. The curves are 
sin^ 6 fits to the data in the range 0:7r/5. (Color Online) 

Fig. [31 showing a deviation from the usual sin^ 9 angular dependence at temperatures close 
to Tc (700 K). In fact a consistent flattening of the free energy is seen in the proximity of the 
hard axis - the high anisotropy energy causes a reduction in the length of the magnetization 
which effectively means it is competing with the exchange interaction. The balance of these 
effects is that the length of magnetization is slightly lower and therefore so is the torque, 
leading to a flattening of the free energy in the maximum. This leads to a lower than 
expected energy barrier which is important for the calculation of relaxation times at high 
temperature. The low temperature scaling of the anisotropic free energy is plotted in Fig. IH 





-0.2 

^ -0.4 

D- 

V- 

2k -0.6 

^ -0.8 













□ 














□ 







-0.5 -0.4 -0.3 -0.2 -0.1 

In (M/Mj) 



FIG. 4: Temperature scaling of anisotropy in LIq FePt. 



showing excellent agreement with the previous theoretical and experimental resulta^i^i^. 
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VI. THIN FILM SYSTEMS WITH NEEL SURFACE ANISOTROPY 



So far we have demonstrated the abihty of the constrained Monte Carlo method to repro- 
duce results which are well known. In the following the method is applied to a system where 
the temperature dependent behavior is unknown - namely thin films with surface anisotropy. 
Understanding surface anisotropy presents a number of challenges due to its complexity, es- 
pecially in nanoparticles^. Due to symmetry thin films present a special case for surface 
anisotropy, where the behavior is much simplified. Thin films have attracted a great deal of 
research interest over the past 50 years and so a large body of experimental data exists^^'^'^. 
Nevertheless, achieving good experimental data on the temperature dependence of surface 
anisotropy requires the creation of very thin films with very sharp interfaces, which has only 
been technologically feasible within the last decade. This is because the influence of surface 
anisotropy is usually determined by varying the thickness of the magnetic layer, so that 
volume and surface contributions can be separated. For thick films the volume component 
strongly dominates the overall anisotropy, leading to a large degree of uncertainty in the 
strength of the surface contribution. Another problem arises with temperature dependent 
atomic migration, structural changes and interface mixing, which cause a change in the 
surface properties'^. 

Calculation of temperature dependent anisotropy 

The constrained Monte Carlo method allows for a thorough investigation of the tem- 
perature dependence of anisotropy in thin films with the Neel surface anisotropy. In the 
case of a perfect single crystal magnetic film with a face-centered-cubic (fee) or simple cubic 
(sc) crystal structure with interfaces cut along the [001] direction, the on-site Neel surface 
anisotropy yields a purely uniaxial anisotropy. In order to simulate a section of such a thin 
film, a generic magnetic material with fee (Tc = ISOO-ft') or sc (Tc = lOOOif) crystal structure 
was chosen. In order to eliminate edge effects within the film, periodic boundary conditions 
in the film plane were used. The surface anisotropy is normally found to be much stronger 
than bulk-type anisotropy, and so a value of = lOKu was chosen. When studying thin 
films with surface anisotropy, a number of basic combinations of anisotropies are possible. 
Principally, in the case of bulk uniaxial anisotropy, the surface and bulk anisotropies can 
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have aligned or opposing easy axes, depending on the sign of the constant. Alternatively, a 
material could possess a cubic bulk anisotropy and uniaxial surface anisotropy. In the fol- 
lowing we present calculations of temperature dependent effects in these thin film systems 
with surface anisotropy. 

Where both bulk and surface anisotropics are uniaxial, the torque curves are similar to 
ones presented in Fig. [H showing a sin(2^) angular dependence. In Fig. |5] we present a 
more complicated situation, where the thin film has cubic bulk anisotropy and Neel surface 
anisotropy. One of the easy axes of the cubic anisotropy coincides with the surface anisotropy 
easy axis (perpendicular to the thin film surface). The torque curve clearly shows a sum- 
mation of uniaxial and cubic anisotropy contributions. The temperature dependence of 
both contributions is presented in Fig. Ei Similar to our results in the previous section, 
the cubic anisotropy exhibits a much stronger temperature dependence than the uniaxial 
(surface) part. Consequently, at low temperature the cubic anisotropy dominates while at 
high temperature the uniaxial surface anisotropy dominates. 
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FIG. 5: The torque curve in a thin film of = 10 atomic planes with fee structure, cubic bulk 
anisotropy and Neel surface anisotropy perpendicular to the thin film plane and parallel to one of 
the bulk anisotropy easy axes for T = lOK. The line provides a guide to the eye. 

Modeling of the reorientation transition in thin films. 

The temperature dependence of the surface anisotropy leads to a number of interesting 
effects, such as a temperature dependent reorientation of the magnetization direction from 
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Temperature [K] 

FIG. 6: Temperature dependence of uniaxial and cubic effective anisotropies in a thin film of 
Lz = 10 atomic planes with fee structure, cubic bulk anisotropy and Neel surface anisotropy 
perpendicular to the thin film plane and parallel to one of the bulk anisotropy easy axes. The lines 
provide a guide to the eye. (Color Online) 

out-of-plane to in-plane and vice versa^""— . Such an effect can occur when the easy directions 
of the surface and bulk anisotropies compete. At low temperatures the magnetization lies 
along the surface easy direction, e.g. perpendicular to the plane. As the temperature 
is increased the surface contribution to the anisotropy energy rapidly decreases, so the 
system magnetization lies along the bulk easy direction, e.g. in the plane. The temperature 
dependence of the effective anisotropy is plotted in Fig. [7] for different thin film thicknesses. 
Given the large difference in the surface and bulk anisotropy constants, the ultrathin films 
fail to show any reorientation transition. As the film thickness is increased the reorientation 
transition becomes more pronounced and occurs at a lower temperature. These results are 
comparable to mean-field calculations by Hucht et alM. 

One other interesting property of the temperature dependent reorientation transition 
is that, depending on the choice of non-magnetic interface material, the temperature of 
the transition can be tuned. To illustrate this phenomenon. Fig. [8] shows a plot of the 
temperature dependence of the total anisotropy for different Neel anisotropy constants, 
emulating the effect of changing the interface material. Here the bulk anisotropy is assumed 
to be perpendicular to the plane, and the surface anisotropy is easy plane. 
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FIG. 7: Temperature dependence of the effective anisotropy in thin films with in-plane bulk 
anisotropy and perpendicular to the plane surface anisotropy for various thin film thicknesses. The 
lines provide a guide to the eye. (Color Online) 
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FIG. 8: Temperature dependence of effective anisotropy for a thin film with competing surface 
(easy plane) and bulk (parallel to z axis) anisotropies for a system of 26 x 26 x 6 unit cells for 
different magnitudes of Kg. The lines provide a guide to the eye. (Color Online) 



Temperature dependence of surface anisotropy 



The temperature dependent effective anisotropy in thin films with surface anisotropy is 
not an intrinsic parameter since it is strongly dependent on the thin film thickness. In 
this section we present two methods which enable the separation of the surface and bulk 
anisotropy contributions as a function of temperature in thin films, for simplicity, with 
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parallel surface and bulk uniaxial anisotropy axes perpendicular to the thin film plane. This 
allows the extraction of the intrinsic uniaxial and surface contributions, independent of the 
thin film thickness. 

The first method is based on the variation of the effective anisotropy with the number of 
surface atoms via the following well known expression^: 

Ko, = Kf + ^{Kf - Kf) (4) 

where Ns and are the number of surface and total atoms, and K^^ and are the 
effective surface and bulk anisotropics, respectively. In Figj9] we present results at different 
temperatures in thin films with a simple cubic lattice and different thicknesses. The data are 
perfectly scaled with the ratio Ng/N and for increasing film thickness the effective anisotropy 
tends towards the temperature dependent bulk value. The fitting of the data to Eq. (jl]) allows 
the extraction of the surface anisotropy constant as a function of temperature, as presented 
in Fig. [ini In this system the surface anisotropy shows a linear decrease with temperature, 
similar to the experimental resulta^i^i^. 

4 
3.5 

3 

^3 2.5 
V 2 
1.5 
1 

0.5 

FIG. 9: Scaling of the effective anisotropy with the ratio between surface and total number of 
atoms Ng/N in a thin film with sc structure, uniaxial bulk anisotropy and Neel surface anisotropy 
both perpendicular to the thin film plane. The lines show the fits to the data. (Color Online) 

Since the surface atoms must be identified in order to calculate the Neel surface 
anisotropy, one can also resolve the surface and bulk contributions to the restoring torque 
curves. Fig. [11] shows the torque curves for the total, bulk and surface parts, evaluated for 
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FIG. 10: Temperature dependence of the effective surface anisotropy in a thin film with sc struc- 
ture, uniaxial bulk anisotropy and Neel surface anisotropy perpendicular to the thin film plane 
determined from the size scaling of the effective anisotropy. The line is a linear fit to the data. 



the same thin film as above. As can be seen, the surface torque also follows the sin (2^) 
behavior and the effective surface anisotropy can therefore be extracted. The values of the 
surface anisotropy obtained through this method are comparable with the ones obtained 
through the scaling method shown in Fig. [TOl 
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FIG. 11: Simulated restoring torque in a thin film of = 10 atomic planes with sc structure, 
uniaxial bulk anisotropy and Neel surface anisotropy perpendicular to the thin film plane for 
T = 10-fC. The lines show a fit of the data to sin 20. (Color Online) 



In practice, the surface torque is a noisy quantity due to the relatively small number of 
atoms. In order to obtain a good thermodynamic average it is generally necessary to use a 
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large number of steps. However, the total torque converges much more rapidly and requires 
relatively few steps. 

The scaling behavior of the anisotropy with magnetization 

The separation of surface and bulk contributions to the anisotropy allows the investigation 
of the temperature scaling of the surface anisotropy separately with respect to the surface 
magnetization, Mgurface, and bulk magnetization, Mbuik- We should bear in mind that the 
magnetization fluctuations on the surface are dependent on the bulk so that in general the 
corresponding scaling exponent is strongly system size dependent. 

For an isolated surface layer the scaling of the surface anisotropy with respect to the sur- 
face magnetization should follow K^^ ~ ^surface' was found in the bulk case. In principle 
this effect could be measured experimentally using a monolayer of magnetic material, though 
such a structure is generally unstable at anything other than cryogenic temperatures. 

We have found that a magnetic thin film with zero bulk anisotropy also follows the 
K^^ ~ Mg^j,fg^^g law, at least for the thin film thicknesses for which our calculations were 
feasible. We have simulated a thin film system with fee crystal structure, surface anisotropy 
perpendicular to the thin film plane and with zero bulk anisotropy. This essentially ensures 
that the only anisotropic contribution to the Hamiltonian comes from the surface. The 
normalized magnetization and surface anisotropy calculated via the torque method as a 
function of temperature are plotted in Fig. [12] for the system dimensions 32 x 32 x 12 unit 
cells. The surface, bulk and volume average magnetization are plotted, each having the 
same Curie temperature but with a different criticality, as previously reported by Binder et 
al^. 

The reduced criticality in the surface magnetization arises from a reduction in coordina- 
tion number. An isolated surface layer would also have a reduced Curie temperature, but in 
our case the surface layer is polarized by the bulk and thus has the same Tc of around 1300 
K. Fig. [T3] shows the temperature scaling of the surface anisotropy with the surface magne- 
tization showing a low temperature exponent of K^^ ~ ^sTirfacc excellent agreement with 
the Callen-Callen theory for single ion uniaxial anisotropy. 

The scaling of the total effective anisotropy in the presence of the Neel surface anisotropy 
with total magnetization is unknown a-priori, and is coordination number, thin film thick- 
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FIG. 12: Plot of normalized magnetization and surface anisotropy against temperature for a thin 
film with zero bulk anisotropy. The surface magnetization shows stronger criticality than the bulk 
and average magnetization. The lines provide a guide to the eye. (Color Online) 
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FIG. 13: Plot of temperature scaling of surface anisotropy with surface magnetization, for a system 
with zero bulk anisotropy. 

ness and material dependent. Nevertheless, it is this scaling which would be measured 
experimentally. To illustrate this effect we have calculated the effective scaling exponent 
K'^ ~ ^2verage ^ systcm with both uniaxial and surface anisotropics for different film 
thicknesses, as shown in Fig. [TH For an isolated surface layer = N the critical exponent 
7 = 3 is recovered. The critical exponent has a maximum and should tend again to 7 = 3 
value for very thick films. 
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FIG. 14: Plot of scaling exponent with thin film thickness in a thin film with Kg = lOK^, parallel 
out-of-plane easy axes and sc lattice. The line provides a guide to the eye. 

VII. CONCLUSION 

We have developed a constrained Monte Carlo method, by which we can compute ther- 
modynamic properties of magnetic systems as a function of the magnetization direction. 
We have shown its novel capability to evaluate the temperature dependence of the magnetic 
anisotropy, which is an important quantity for technological applications in hard magnetic 
materials. The method has been utilized to compute the temperature dependence of bulk 
magnetic anisotropy and we have recovered numerically the analytic scaling law of Callen 
and Callen. 

The importance of the method resides in its potential to calculate temperature-dependent 
effective anisotropics in complex materials. The present challenge for modeling of magnetic 
materials is the multiscale approach, where the ab-initio information is passed to a different 
scale with the aim to model larger sizes of material, using, for example micromagnetics. 
The CMC method provides a real possibility to link the quantum mechanical scale with 
micromagnetics, via the parameterization of a classical Heisenberg Hamiltonian^^. This 
Hamiltonian can be used to calculate temperature dependent equilibrium behavior. To 
demonstrate this, we have applied the method to the calculation of the effective anisotropy 
bulk FePt, using the model parameterized through ab-initio calculations, and have recovered 
the experimentally observed -f^FcPt ~ M^'^ temperature scaling. Moreover, we have also 
noted the appearance of a flattening of the free energy surface at the energy maximum, an 
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effect not previously seen with less subtle methods. 

We have then applied the method to a variety of thin films with surface anisotropy, inves- 
tigating the size dependent behavior and temperature dependence of surface anisotropy. We 
have shown the capability of the method to simulate the temperature-dependent magneti- 
zation reorientation transition. We also have shown that the method enables the separation 
of the temperature dependent surface anisotropy as an intrinsic system parameter, indepen- 
dent of the thin film thickness. Our results demonstrate a linear temperature dependence 
of the surface anisotropy, consistent with the experimental results in Gdr^^, Ni^ and Fe^^ 
grown on different substrates. However, when comparing with our results, various factors 
should be taken into account, including structural changes with increased temperatures or a 
lattice mismatch which could influence the bulk anisotropy. One other possibility is that of 
an enhanced exchange interaction at the surface of the material^. An increased exchange 
interaction would lead to an increase in the criticality of the surface layer and commensu- 
rate reduction in the temperature dependence of the surface anisotropy. In the future, more 
complicated models could take into account these effects. 

Finally, we have investigated the scaling of the anisotropy in thin films with magnetiza- 
tion. In thin films with zero bulk anisotropy, the surface anisotropy scales with the surface 
magnetization following the Callen-Callen law. In other cases we report no universal scaling 
behavior of the surface anisotropy. 

To summarize, the constrained Monte Carlo method is a powerful tool, allowing to include 
both thermodynamic fluctuations and entropy into the evaluation of macroscopic quantities 
such as temperature dependent magnetic anisotropy. In the future we plan to apply the 
method to model more realistic systems within a general trend to large-scale material mod- 
eling aiming the design of novel materials with potential applications. 
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Appendix A: Jacobians in the trial moves 



We are displacing two spins, whose combined volume element is d'^Si(PSj. However, we 
plan to readjust Sj to keep M constant. Accordingly, we eliminate the variable Sj in favor of 
M and rewrite the volume element as | J(Sj, M)| d'^Si d'^M, expressing Sj in terms of M and 
expecting the Jacobian J {Si, M) to be nontrivial. This allows us to view the trial move as 
taking place in the Sj, M variables. The Boltzmann probability density in these variables is 
p oc I J| exp{—f3'H). Since our trial moves are reversible in (Sj, M), the Metropolis-Hastings 
test ratio^ is just p'/ p = \J' / J\ exp(— /^AH) . The fact that we always keep M' = M has no 
bearing on the argument and we use the same ratio to decide whether to accept the move 
from (Si,M) to (S^,M). 

The first result we need is the relationship between a spherical surface element and its 
projection in the XY plane. We have 

_ dSjx dSjy d^S- = ^^^^ ^^jy (Al) 

* I Q I ' I Q I 

d'M=^^^^^ = dM,dMy . (A2) 
|M,| 

Too see this, note that for any unit vector S the angle between the tangent plane to the 
sphere at S and the XY plane is simply 9 = cos~^ (^I'S'zlj, therefore the projected area 
dSxdSy is |cos(6')| d'^S = \Sz\ d^S. Equations f lAip and flA2p are the inverse of this relation, 
stated for the three vectors Sj, Sj, M. In the case of M, Eq. (lA2p . the denominator \Mz\ is 
not necessary because we are locking M in the positive z direction, therefore Mz = 1. 

Next, we express dM^ and dMy in terms of dSix, dSiy, dSjx and dSjy. We assume that 
M lies in the z direction but we do not make the same assumption for M + dM.. For the 
nnnormalized vector M, we have 

dMx = dSix + dSjx , dMy = dSiy + dSjy (A3) 

since the other spins are fixed. For the normalized M, we have 

d'M = d (M-^M) = dm - iVf-^M dM . (A4) 



We only need the x and y component of flA4p . Since M lies in the +z direction, the last 
term contributes nothing to the in-plane components and we don't need to calculate dM 
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(although it would be easy to do so). The remaining term gives 



dMx = M-^ dMx , dMy = M"^ dMy 



(A5) 



We substitute (jM]) in (^MS> and replace M by M^: 

dMx = M;^ dSix + M-^ dSjx 



dMy = m:^ dSiy + m:^ ds 



(A6) 



We rewrite flA6P in matrix form, 



dSix 




1 













dSix 


dSiy 







1 










dSiy 


dMx 
















dSjx 


dMy 
















_dSjy 



(AT) 



The Jacobian of the change of variables is the determinant of the matrix in flA7l) . namely 
M~^. Therefore, the volume elements are related by 



dSix dSiy dli/lj^ d]\d^y — z ^^ix dSiy dSjx ^^jy 
and the inverse relation is 



dSix dSiy dSjx dSjy = M^ dSix dSiy dMx dMy 



Dividing by \Siz \ \Sjz\ and using (lAip . we obtain 

d Si d S j 



2^ -^20 _ £s,d'M 



\S. 



(A8) 



(A9) 



(AlO) 



Thus the Jacobian is \J\ = M'^/\Sjz\ and we accept each trial move with probability 



P = min [1, \J'/J\ exp(-/3AH)] 



mm 



l,l^)giexpMA«) 



(All) 



as stated in step [7] of the algorithm. 
Appendix B: Ergodicity 



We show that every admissible state (i.e. with total magnetization along +z) can be 
reached from any other admissible state by a sequence of constrained Monte Carlo trial 
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moves. In this context we may ignore the Metropohs acceptance test at step [8] of the 
algorithm, but the trial moves must still pass the kinematic tests at steps H] and [51 The 
proof relies on two lemmas: 

Lemma 1 Any spin with a negative z component can be moved to the hemisphere Sz > 
by a sequence of trial moves, without any other spin crossing the plane z = 0. 

Lemma 2 Assuming all the spins are in the hemisphere Sz > 0, one of them can be moved 
to the z axis by a sequence of trial moves, without any spin crossing the plane z = 0. 

Repeated applications of Lemma [T] allow us to move the spins one by one to the hemi- 
sphere Sz > 0, at which point Lemma [2] becomes apphcable. 

Repeated applications of Lemma |5] allow us to move all the spins to the positive z axis. 
After the first application we have one spin along +z; the remaining A^— 1 spins are perturbed 
in the process, but they remain in the upper hemisphere and their final magnetization must 
also point along +z. That is, they form an admissible state of A^ — 1 spins. Lemma |2]becomes 
applicable to them and we can iterate. Thus, any admissible state can be collapsed to the 
saturated state by a sequence of trial moves. By chaining such a collapsing sequence with 
the inverse of another one, we can connect any two admissible states. 

Proof of Lemma [2] 

We use trial moves where the primary spin does not cross the plane z = 0. We write Sj 
(lower case) for the projection {Si^, Siy) in the XY plane. It is enough to consider the Sj 
since the z components are uniquely determined by them, Siz = +a/1 — ||sjP. The trial 
moves can not fail at step El leaving only step H] to consider. A trial move then consists in 
displacing two points Sj, s^- by equal and opposite amounts, while keeping both points in the 
unit disk, ||sj||, ||Sj|| < 1. 

We want to show that, for any set of A^ points in the unit disk with centroid (0, 0), one of 
the points can be shifted to (0, 0) by a combination of such moves. In fact we have a slightly 
more general result: in any set of A^ points, one can be moved to the centroid, whether or 
not the centroid is (0,0). 

We proceed by induction on A^. For N < 2 there is nothing to prove. For N = 2, let 
c = (si + S2)/2 be the centroid. Since c is invariant, the condition ||s2|| < 1 can be expressed 
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in terms of Si, ||2c — Si|| < 1. Similarly we have ||2c — S2II < 1. This means that we can move 
the two points symmetrically about c, as long as we keep them both within a lens-shaped 
"admissible zone", which is the intersection of the unit disk and of its inversion through c. 
This is shown in Fig. [T3 The admissible zone, being convex, contains c. Therefore we can 
move both points to c, using several sub-moves if necessary. 



FIG. 15: Admissible zone for projected trial moves in the proof of Lemma [2j The projections of 
the two spins (black circles) can move symmetrically about their centroid (white circle) within the 
region shown in white. 

For > 2 we assume by the induction hypothesis that the centroid Cjv-i of the first A^ — 1 
points is already occupied by one of the s^. The full centroid c^v = (s^r + {N — l)sfc)/A^ 
lies on the segment [sfc,SAr] and is in the admissible zone of Sk,SN- Therefore we can move 
either point to the centroid. This proves the lemma. 

Proof of Lemma [1] 

We want to reduce by one the number of spins with 5*2 < 0. Our rules allow trial moves 
where the primary spin flips its z component and the secondary spin does nothing, but we 
would like to avoid these potentially low probability moves. Our strategy is to move the 
offending spin close to the plane 2; = so as to require only a small jump in Sz in the final 
move. We proceed as in the proof of Lemma [21 using projections in the XY plane, but this 
time we have to avoid rejection at step [5] of the algorithm. 
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Consider an admissible state where spin i has Si^ < 0. Choose a j such that Sj^ > 
(there must be one, since is positive). We use trial moves with Sj as the primary spin 
and Sj as the compensation spin. Let Sj, Sj be the projections in the XY plane. The two 
points Sj, Sj are restricted to a lens-shaped admissible region as in the proof of Lemma [2l 
but there may be a further restriction to ensure that remains positive. 

As before, we express Sj in terms of Sj, Sj = 2c — Sj. The contribution of spins i and j to 
Mz is then 

Si, + Sj, = - VT^NF + Vl- ||2c-s,||2 . (Bl) 

It can be shown that the contours of (IBip in the {six, Siy) plane are arcs of ellipses centered 
on c and tangent to the admissible zone, as shown in Fig. [161 and that f IBip increases 
monotonically as Sj moves from the interior edge to the exterior edge of the admissible zone. 
If Sj stays on the distal side of its starting contour, the value of cannot become negative. 




FIG. 16: Admissible zone for projected trial moves in the proof of Lemma [TJ The dashed lines are 
contom's of Si, + Sj, as a function of Six, Siy. A possible path for the projection of the primary 
spin is shown in black, with the matching path of the secondary spin in grey. If the primary spin 
stays in the white region the total magnetization M, cannot become negative. 

This is represented by the white region in Fig. [161 We see that there is an ample supply of 
paths that take Sj arbitrarily close to the edge of the unit disk. At that point a final trial 
move can take Sj across the plane Si, = 0. The only other spin involved, Sj, remains in the 
positive z hemisphere throughout. This proves Lemma [H 



26 



Appendix C: Macroscopic Torques 



The magnetization direction M plays the role of a macroscopic parameter. In this section 
we identify the generalized forces acting on M (torques) with the derivatives of the Helmholtz 
free energy. 

By way of introduction, consider a thermodynamic system with an external macroscopic 
parameter q and microscopic states S,a- We treat a as a discrete index for simplicity. The 
energy of each microscopic state is a function of g; the partition function and the Helmholtz 
free energy are — 

Z{q) = Y,^^Y>{-m{Uq)) , (CI) 

a 

^(g) = -/3-^Ln(Z(g)) . (C2) 
Differentiation of (lC2p with respect to q yields immediately 

a 

= {-m/dq) 



(C3) 



where in the second line of flC3P we recognized the sum over microstates as a thermodynamic 
average, since the weights exp(— /^^/(^q, g)) are the Boltzmann probabilities. Thus the 
mean force conjugate to q is the negative derivative with respect to q of the Helmholtz free 
energy F{q). 

The right-hand side of ( ]C3p . being a thermodynamic average, can be computed by the 
Metropolis algorithm. As a general rule, derivatives and differences of free energies are 
computable in this way, even though the free energy itself is not. 

The argument is not directly applicable to our case since M is not a parameter in the 
Hamiltonian, but a restriction on the set of admissible states. One could in principle set 
up a system of 2N — 2 coordinates that forms a complete set with M and treat M as a 
parameter, but we never did this. Instead, we picked a temporary coordinate system for 
each Monte Carlo move, to be discarded after the move. Indeed, the whole of appendix |X] 
is a stratagem to avoid setting up global coordinates. 

However we can recover a result similar to (1C3I) . Given any two directions M and M' 
we can find a rotation TZ that sends M to M'. If we apply this rotation globally to every 
spin Sj that constitutes a microstate we obtain a measure-preserving bijection between 
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the two admissible manifolds. This allows us to replace sums over the microstates of the 
M' manifold by sums over the microstates (rotated) of the M manifold. In particular, the 
partition function Z' = 2(M') is 

Z' = 5^exp(-/3H(7^U) • (C4) 

a 

Dividing by Z and factoring out an exp(— /3'H(.^q,)) from each term, 
Z7Z = ^Z-iexp(-/3-H(U)x 

a 

exp[-/3{'Hin^a)-'H{U)] , (C5) 



or 



z'/z = (^exp[-(]{'H{no-n{0)]) , (ce) 

where the sums in (1C4HC5P and the thermodynamic average in (]C6p are over the admissible 
manifold of M. Taking logarithms and dividing by — /3, we have 

r-T = -r' Ln [(exp [-p{n{no - m))] )] (C7) 

where again the thermodynamic average is over the manifold of M. 

Specializing to an infinitesimal rotation, TZv = v + dOh x \ + 0{d6'^) and writing each mi- 
crostate as a collection of spins Sj, the energy difference in the argument of the exponential 
is 

H(7^0 - m) = rf^^ $^(n X S.) • -4^ + 0{d9') (C8) 



deh■J2(^^ X ^) +0id9^) . 
■ \ dSi / 



(C9) 



As d6 approches zero the logarithm and exponential in ( \C7\i become a no-op, 

jr'_jr= rf^n. ^^Si X +0(^0') . (CIO) 

Meanwhile, M' is related to M by the same infinitesimal rotation and we can expand J-"' = 
J-'(M') in a Taylor series. 

= T{M + dehxM + 0{de'^)) (Cll) 

= T + d9(hxM) ■ ^ + 0(d9^) (C12) 

= + dOh ■ (m X + 0{de^) (C13) 
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Combining the terms of order d9 in ( ]C10|) and (IClSp . we obtain 



X 



m 



n ■ M X 



dm) 



(CM) 



Since the equahty holds for all fi, we have finally 

dU 



X 



M X 



dm 



(C15) 
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